Differential abundance and differential expression analysis in cardiomyocytes

workflowr

• Summary
• Checks
• Past versions

Last updated: 2024-08-20

Checks: 7 0

Knit directory: DCM_snRNAseq/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20240606)

The command set.seed(20240606) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: c979c9e

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version c979c9e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/figure/

Untracked files:
    Untracked:  analysis/omnipathr-log/
    Untracked:  code/Add_metadata.R
    Untracked:  code/Check DYSF.R
    Untracked:  code/Clustering_genes.R
    Untracked:  code/DE_5_percent.Rmd
    Untracked:  code/DE_5_percent.html
    Untracked:  code/DE_5_percent/
    Untracked:  code/DE_CM_test1.R
    Untracked:  code/DE_no_401.Rmd
    Untracked:  code/DE_no_401.html
    Untracked:  code/DE_no_401/
    Untracked:  code/Differential abundance test1.R
    Untracked:  code/Differential_Expression_edgeR_All.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_2.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_2.html
    Untracked:  code/Differential_Expression_edgeR_All_Age.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_Age_2.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_Age_2.html
    Untracked:  code/Differential_Expression_edgeR_All_groups.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_groups.html
    Untracked:  code/Differential_Expression_edgeR_All_groups_2.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_groups_2.html
    Untracked:  code/Differential_Expression_edgeR_All_groups_3.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_groups_3.html
    Untracked:  code/PCA and CCA analysis test.R
    Untracked:  code/Pseudobulk DCM HC.R
    Untracked:  code/QC_integration_annotation_orig.Rmd
    Untracked:  code/UpSet_plot_DEGs.Rmd
    Untracked:  code/Volcano_highlighted_genes.R
    Untracked:  code/ezInteractiveTable.Rmd
    Untracked:  code/ezInteractiveTable.html
    Untracked:  code/old.R
    Untracked:  core
    Untracked:  data/Cellbender_output/
    Untracked:  data/Cellranger_output/
    Untracked:  data/DCM_Clinical_data.xlsx
    Untracked:  data/DCM_Clinical_data_26.xlsx
    Untracked:  data/Raw/
    Untracked:  omnipathr-log/
    Untracked:  output/Cardiomyocytes_DA_DE/
    Untracked:  output/Cardiomyocytes_DA_DE_26/
    Untracked:  output/Cardiomyocytes_subclustering/
    Untracked:  output/Cardiomyocytes_subclustering_26/
    Untracked:  output/Differential_expression_edgeR_All/
    Untracked:  output/Differential_expression_edgeR_All_Age/
    Untracked:  output/QC_integration_annotation/
    Untracked:  output/QC_integration_annotation_26/

Unstaged changes:
    Modified:   analysis/Cardiomyocytes_DA_DE.Rmd
    Deleted:    analysis/Differential_Expression_edgeR_All.Rmd
    Deleted:    analysis/Differential_Expression_edgeR_All_Age.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Cardiomyocytes_DA_DE_26.Rmd) and HTML (docs/Cardiomyocytes_DA_DE_26.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	c979c9e	GinoBonazza	2024-08-20	wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd")
html	182f9de	GinoBonazza	2024-08-19	Build site.
html	5cd0f74	GinoBonazza	2024-08-01	Build site.
Rmd	f3f2075	GinoBonazza	2024-08-01	wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd")
html	61e8d17	GinoBonazza	2024-08-01	Build site.
Rmd	9e4b482	GinoBonazza	2024-08-01	wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd")

---

Setup

# Get current file name to make folder
current_file <- "Cardiomyocytes_DA_DE_26"

# Load libraries
library(here)
library(readr)
library(readxl)
library(xlsx)
library(Seurat)
library(DropletUtils)
library(Matrix)
library(scDblFinder)
library(scCustomize)
library(dplyr)
library(ggplot2)
library(magrittr)
library(tidyverse)
library(reshape2)
library(S4Vectors)
library(SingleCellExperiment)
library(pheatmap)
library(png)
library(gridExtra)
library(knitr)
library(scales)
library(RColorBrewer)
library(Matrix.utils)
library(tibble)
library(ggplot2)
library(scater)
library(patchwork)
library(statmod)
library(ArchR)
library(clustree)
library(harmony)
library(gprofiler2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(AnnotationHub)
library(ReactomePA)
library(statmod)
library(edgeR)
library(speckle)
library(EnhancedVolcano)
library(decoupleR)
library(OmnipathR)
library(dorothea)
library(enrichplot)
library(png)
library(ezRun)
library(UpSetR)
library(ComplexHeatmap)

#Output paths
output_dir_data <- here::here("output", current_file)
if (!dir.exists(output_dir_data)) dir.create(output_dir_data)

if (!dir.exists(here::here("docs", "figure"))) dir.create(here::here("docs", "figure"))

output_dir_figs <- here::here("docs", "figure", paste0(current_file, ".Rmd"))
if (!dir.exists(output_dir_figs)) dir.create(output_dir_figs)

Load cardiomyocytes dataset

CM <- readRDS(here::here("output", "Cardiomyocytes_subclustering_26", "CM.rds"))

DefaultAssay(CM) <- "RNA"
CM <- NormalizeData(CM)

Differential abundance analysis

#Prepare matadata for differential abundance and differential expression testing
metadata <- CM@meta.data %>%
  dplyr::select(Sample:CO_l.min) %>%
  unique() %>%
  dplyr::select(-LVID_mm, -RVIT_mm) 
metadata <- metadata %>%
  dplyr::select(1:3, 5, 6, 4, 7:ncol(metadata)) %>%
  dplyr::arrange(match(Sample, names(table(CM$Sample))))

rownames(metadata) <- metadata$Sample

# Check if metadata$Sample matches names(table(CM$Sample))
all.equal(metadata$Sample, names(table(CM$Sample)))

[1] TRUE

ezInteractiveTableRmd(metadata)

# Filter for Group == "DCM"
metadata_DCM <- metadata %>%
  dplyr::filter(Group == "DCM")

# Initialize an empty data frame for storing differential abundance results
differential_abundance <- data.frame()
props <- getTransformedProps(CM$cell_state, CM$Sample, transform="logit")

# Loop through each metadata column from Age_years to CO_l.min
for (i in 6:(length(metadata_DCM))) {
  metadata_subset <- metadata_DCM[complete.cases(metadata_DCM[[i]]),]
  props_subset <- props
  props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
  props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
  props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
  metadata_subset <- metadata_subset %>% 
    arrange(Sample, colnames(props_subset$Proportions))
  design <- model.matrix(~ metadata_subset[[i]])
  fit <- lmFit(props_subset$TransformedProps, design)
  fit <- eBayes(fit, robust=TRUE)
  differential_abundance_temp <- topTable(fit, n = Inf, coef = 2)
  differential_abundance_temp$metadata <- colnames(metadata_subset[i])
  differential_abundance_temp$cell_state <- rownames(differential_abundance_temp)
  differential_abundance <- rbind(differential_abundance, differential_abundance_temp)
}

# Clean up the environment
rm(metadata_subset, props_subset, design, fit, differential_abundance_temp)

# Reset row names and arrange the results by adjusted p-value
rownames(differential_abundance) <- NULL
differential_abundance <- dplyr::arrange(differential_abundance, adj.P.Val)

write.csv(differential_abundance, file = here::here(output_dir_data, "CM_Differential_Abundance.csv"), quote=F, row.names = F)

t <- ezInteractiveTableRmd(differential_abundance)
t[["x"]][["options"]][["pageLength"]] <- 10
t

metadata_subset <- metadata[complete.cases(metadata[["LVEDD_mm"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$LVEDD_mm)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$LVEDD_mm, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="LVEDD_mm", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Past versions of Speckle_LVEDD_mm-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_subset, props_subset, design, fit.prop)

metadata_subset <- metadata[complete.cases(metadata[["LVEF_percent"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$LVEF_percent)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$LVEF_percent, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="LVEF_percent", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Past versions of Speckle_LVEF_percent-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_subset, props_subset, design, fit.prop)

metadata_subset <- metadata[complete.cases(metadata[["Age_years"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$Age_years)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$Age_years, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="Age_years", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Past versions of Speckle_Age_years-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_subset, props_subset, design, fit.prop)

metadata_subset <- metadata[complete.cases(metadata[["RVSP_mmHg"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$RVSP_mmHg)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$RVSP_mmHg, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="RVSP_mmHg", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Past versions of Speckle_RVSP_mmHg-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_subset, props_subset, design, fit.prop)

design <- model.matrix(~ 0 + metadata$Group)
colnames(design) <- c("DCM", "HC")
mycontr <- makeContrasts(DCM-HC, levels=design)
differential_abundance_HC <- propeller.ttest(props, design, contrasts = mycontr, robust=TRUE, trend=FALSE, sort=TRUE)

differential_abundance_HC

    PropMean.DCM PropMean.HC PropRatio Tstatistic      P.Value         FDR
CM3   0.11135354  0.03836732 2.9023016   3.800582 0.0006739595 0.003369797
CM4   0.08662428  0.02831785 3.0589987   3.202724 0.0032613900 0.008153475
CM1   0.50911696  0.57790022 0.8809773  -1.527054 0.1374288591 0.229048098
CM2   0.24032230  0.29095689 0.8259722  -1.002228 0.3244063073 0.405507884
CM5   0.05258291  0.06445772 0.8157738  -0.479165 0.6353640356 0.635364036

par(mfrow = c(1, 5), mar = c(3, 5, 3, 1))
for (i in seq(1, 5, 1)) {
  stripchart(props$Proportions[i,] ~ metadata$Group,
             vertical = TRUE, pch = 16, method = "jitter", 
             col = c("orange", "purple"),
             cex = 2, ylab = "Proportions", cex.lab = 2, cex.axis = 2)
  title(rownames(props$Proportions)[i], cex.main = 2.5)
}

Past versions of Speckle_Group-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(design, mycontr)

Differential expression analysis

Remove genes expressed in <1% of cells and MT and RP genes

percent_stats <- Percent_Expressing(seurat_object = CM, features = rownames(CM), assay = "RNA", entire_object = TRUE)
percent_stats <- percent_stats %>%
  dplyr::rename(percentCells = All_Cells)
percent_stats$gene <- rownames(percent_stats)

write.csv(percent_stats, file = here::here(output_dir_data, "CM_percent_stats.csv"), quote=F, row.names = T)

percent_stats <- read.csv(here::here(output_dir_data, "CM_percent_stats.csv"))
rownames(percent_stats) <- percent_stats$gene

keep_genes <- rownames(dplyr::filter(percent_stats, percentCells > 1))
data <- CM[which(rownames(CM) %in% keep_genes),]
rm(keep_genes)

t <- ezInteractiveTableRmd(percent_stats)
t[["x"]][["options"]][["pageLength"]] <- 10
t

Create and prepare pseudobulk object for differential expression analysis

pseudocounts <- Seurat2PB(data, sample="Sample", cluster = "cell_type")
colnames(pseudocounts) <- pseudocounts$samples$sample

saveRDS(pseudocounts, 
        here::here(output_dir_data, "CM_pseudocounts.rds"))

keep_samples <- pseudocounts$samples$lib.size > 5e4
pseudocounts <- pseudocounts[, keep_samples]
keep_genes <- filterByExpr(pseudocounts)
pseudocounts <- pseudocounts[keep_genes, ]
rm(keep_samples, keep_genes, data)

Create empty objects for the DE output

results <- list()
signif <- list()
volcano <- list()
n_de_genes <- data.frame()

pseudocounts_DCM <- pseudocounts[, metadata_DCM$Sample]
all.equal(metadata_DCM$Sample, colnames(pseudocounts_DCM$counts))

[1] TRUE

for (i in 6:(length(metadata_DCM))) {
  metadata_subset <- metadata_DCM[complete.cases(metadata_DCM[[i]]),]
  metadata_subset[,i] <- rescale(metadata_subset[,i])
  design <- model.matrix(~ metadata_subset[[i]])
  colnames(design)[2] <- colnames(metadata_subset[i])
  
  pseudocounts_subset <- pseudocounts_DCM[, metadata_subset$Sample]
  pseudocounts_subset <- normLibSizes(pseudocounts_subset)
  pseudocounts_subset <- estimateDisp(pseudocounts_subset, design)
  
  fit <- glmQLFit(pseudocounts_subset, design, robust=TRUE)
  fit <- glmQLFTest(fit, coef = 2)
  
  results[[i-5]] <- topTags(fit, n = Inf)$table
  results[[i-5]] <- merge(results[[i-5]], percent_stats[, c("gene", "percentCells")], by = "gene", all.x = FALSE)
  signif[[i-5]] <- results[[i-5]] %>% dplyr::filter(FDR < 0.05, abs(logFC) > 0.5) %>% 
    dplyr::arrange(FDR)

  names(results)[i-5] <- colnames(metadata_subset[i])
  names(signif)[i-5] <- colnames(metadata_subset[i])
  
  write.csv(results[[i-5]], file = here::here(output_dir_data, paste0("CM_DE_Results_", colnames(metadata[i]), ".csv")), quote=F, row.names = F)
  write.csv(signif[[i-5]], file = here::here(output_dir_data, paste0("CM_DE_Significant_", colnames(metadata[i]), ".csv")), quote=F, row.names = F)
  
  n_de_genes_temp <- data.frame(metadata = colnames(metadata_subset[i]),
                                n_upregulated = sum(signif[[i-5]]$logFC > 0.5),
                                n_downregulated = sum(signif[[i-5]]$logFC < -0.5),
                                n_tot = nrow(signif[[i-5]])
                                )
  n_de_genes <- rbind(n_de_genes, n_de_genes_temp)
  
  volcano[[i-5]] <- EnhancedVolcano(results[[i-5]],
                  lab = results[[i-5]]$gene,
                  x = "logFC",
                  y = "FDR",
                  labSize = 0,
                  titleLabSize = 16,
                  subtitleLabSize = 14,
                  axisLabSize = 12,
                  captionLabSize = 9,
                  pointSize = 0.5,
                  FCcutoff = 0.5,
                  pCutoff  = 0.05,
                  ylim = c(0, 3.5),
                  col = c("black", "black", "black", "red"),
                  colAlpha = 1,
                  drawConnectors = FALSE,
                  subtitle = NULL,
                  title = paste0("Cardiomyocytes", "\n", colnames(metadata_subset[i]), "\nn = ", nrow(metadata_subset))
  ) + theme(legend.position = "none")
  names(volcano)[i-5] <- colnames(metadata_subset[i])
}
rm(metadata_subset, design, pseudocounts_subset, fit, n_de_genes_temp)

print((volcano[[1]] | volcano[[2]] | volcano[[3]] | volcano[[4]] | volcano[[5]] | volcano[[6]] | volcano[[7]]) / 
  (volcano[[8]] | volcano[[9]] | volcano[[10]] |  volcano[[11]] | volcano[[12]] | volcano[[13]] | volcano[[14]]))

Past versions of Volcano_all_metadata-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

gene_sets <- list()
for (i in 1:length(signif)) {
  gene_sets[[i]] <- signif[[i]]$gene
  names(gene_sets)[i] <- names(signif)[i]
  }
m <- make_comb_mat(gene_sets, mode = "intersect", min_set_size = 10)
m <- m[comb_degree(m) >= 2]
m <- m[comb_size(m) >= 10]
upset_plot <- UpSet(m, column_title = "Cardiomyocytes", comb_order = order(comb_size(m)), 
                    top_annotation = upset_top_annotation(m, add_numbers = TRUE),
                    right_annotation = upset_right_annotation(m, add_numbers = TRUE))
draw(upset_plot)

Over-representation analysis

# Calculate gene counts across cells
gene_counts <- rowSums(CM@assays$RNA@counts > 0)

# Filter genes to include only those expressed in at least 3 cells
universe_genes <- names(gene_counts[gene_counts >= 3])
write.csv(as.data.frame(universe_genes), here::here(output_dir_data, "CM_universe_genes.csv"), quote=F, row.names = F)

# Convert universe gene symbols to Entrez IDs
universe_entrez <- bitr(universe_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
write.csv(as.data.frame(universe_entrez), here::here(output_dir_data, "CM_universe_entrez.csv"), quote=F, row.names = F)

rm(gene_counts)

up_genes <- filter(signif[["mPAP_mmHg"]], logFC > 0.5)$gene
up_genes <- bitr(up_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
down_genes <- filter(signif[["mPAP_mmHg"]], logFC < -0.5)$gene
down_genes <- bitr(down_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID

reorder_GO_by_pvalue <- function(enrichGO_result) {
  go_results <- enrichGO_result@result
  go_results_sorted <- go_results[order(go_results$p.adjust), ]
  enrichGO_result@result <- go_results_sorted
  return(enrichGO_result)
}

mPAP_GO_up <- enrichGO(up_genes, OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
mPAP_GO_up <- gsfilter(mPAP_GO_up, by = 'Count', min = 3)
mPAP_GO_up <- reorder_GO_by_pvalue(mPAP_GO_up)
saveRDS(mPAP_GO_up, here::here(output_dir_data, "CM_ORA_GO_mPAP_up.rds"))

mPAP_GO_down <- enrichGO(down_genes, OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
mPAP_GO_down <- gsfilter(mPAP_GO_down, by = 'Count', min = 3)
mPAP_GO_down <- reorder_GO_by_pvalue(mPAP_GO_down)
saveRDS(mPAP_GO_down, here::here(output_dir_data, "CM_ORA_GO_mPAP_down.rds"))

mPAP_GO_up <- readRDS(here::here(output_dir_data, "CM_ORA_GO_mPAP_up.rds"))
mPAP_GO_down <- readRDS(here::here(output_dir_data, "CM_ORA_GO_mPAP_down.rds"))
p1 <- dotplot(mPAP_GO_up, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
p2 <- dotplot(mPAP_GO_down, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nLow mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1|p2)

mPAP_GO_simplified_up <- simplify(
  mPAP_GO_up,
  cutoff = 0.75,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
mPAP_GO_simplified_down <- simplify(
  mPAP_GO_down,
  cutoff = 0.75,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
p1 <- dotplot(mPAP_GO_simplified_up, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
p2 <- dotplot(mPAP_GO_simplified_down, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nLow mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1 | p2)

mPAP_KEGG_up <- enrichKEGG(up_genes, organism = 'hsa', universe = universe_entrez)
mPAP_KEGG_up <- gsfilter(mPAP_KEGG_up, by = 'Count', min = 3)
saveRDS(mPAP_KEGG_up, here::here(output_dir_data, "CM_ORA_KEGG_mPAP_up.rds"))
mPAP_KEGG_down <- enrichKEGG(down_genes,  organism = 'hsa', universe = universe_entrez)
mPAP_KEGG_down <- gsfilter(mPAP_KEGG_down, by = 'Count', min = 3)

mPAP_KEGG_up <- readRDS(here::here(output_dir_data, "CM_ORA_KEGG_mPAP_up.rds"))
p1 <- dotplot(mPAP_KEGG_up, showCategory = 10, title = paste0("KEGG Pathways", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1)

mPAP_REACTOME_up <- enrichPathway(up_genes, organism = 'human', universe = universe_entrez, readable = TRUE, pvalueCutoff = 0.1)
mPAP_REACTOME_up <- gsfilter(mPAP_REACTOME_up, by = 'Count', min = 3)
saveRDS(mPAP_REACTOME_up, here::here(output_dir_data, "CM_ORA_REACTOME_mPAP_up.rds"))
mPAP_REACTOME_down <- enrichPathway(down_genes, organism = 'human', universe = universe_entrez, readable = TRUE, pvalueCutoff = 0.1)
mPAP_REACTOME_down <- gsfilter(mPAP_REACTOME_down, by = 'Count', min = 3)

mPAP_REACTOME_up <- readRDS(here::here(output_dir_data, "CM_ORA_REACTOME_mPAP_up.rds"))
p1 <- dotplot(mPAP_REACTOME_up, showCategory = 10, title = paste0("REACTOME", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1)

rm(mPAP_GO_up, mPAP_GO_down, mPAP_GO_simplified_up, mPAP_GO_simplified_down, mPAP_KEGG_up, mPAP_KEGG_down, mPAP_REACTOME_up, mPAP_REACTOME_down)

Gene set enrichment analysis

ranked_genes <- results[["mPAP_mmHg"]]
ranked_genes$metric <- ranked_genes$logFC*-log10(ranked_genes$PValue)
entrezid <- bitr(ranked_genes$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db", drop = TRUE)
ranked_genes <- merge(ranked_genes, entrezid, by.x = "gene", by.y = "SYMBOL")
ranked_genes <- ranked_genes[!duplicated(ranked_genes$metric), ]
genelist <- ranked_genes$metric
names(genelist) <- ranked_genes$ENTREZID
genelist <- genelist[order(genelist, decreasing = T)]
write.csv(ranked_genes, file = here::here(output_dir_data, "CM_ranked_genes.csv"), quote=F, row.names = F)
t <- ezInteractiveTableRmd(ranked_genes)
t[["x"]][["options"]][["pageLength"]] <- 10
t

gseaGO_all <- gseGO(geneList = genelist,
                OrgDb = org.Hs.eg.db,
                ont = "ALL",
                keyType = "ENTREZID",
                nPermSimple = 1000000,
                minGSSize = 10)
gseaGO_all <- reorder_GO_by_pvalue(gseaGO_all)

saveRDS(gseaGO_all, here::here(output_dir_data, "CM_GSEA_GO_mPAP.rds"))

gseaGO_all <- readRDS(here::here(output_dir_data, "CM_GSEA_GO_mPAP.rds"))
gseaGO_all_simplified <- simplify(
  gseaGO_all,
  cutoff = 0.75,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
p1 <- dotplot(gseaGO_all_simplified, showCategory = 10, split=".sign", title = paste0("GO Enrichment Analysis", "\n mPAP-associated genes"), label_format = 30, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"),
        strip.text = element_text(face = "bold", size = 18)) +
  facet_grid(.~.sign)
print(p1)

classify_terms <- function(df) {
  df$sign <- ifelse(df$NES > 0, "activated", "suppressed")
  return(df)
}
gseaGO_all_simplified@result <- classify_terms(gseaGO_all_simplified@result)
# Split the GSEA results
activated_terms <- subset(gseaGO_all_simplified@result, sign == "activated")
suppressed_terms <- subset(gseaGO_all_simplified@result, sign == "suppressed")

# Create two separate enrichResult objects
gseaGO_all_activated <- gseaGO_all_simplified
gseaGO_all_activated@result <- activated_terms
gseaGO_all_activated <- pairwise_termsim(gseaGO_all_activated)

gseaGO_all_suppressed <- gseaGO_all_simplified
gseaGO_all_suppressed@result <- suppressed_terms
gseaGO_all_suppressed <- pairwise_termsim(gseaGO_all_suppressed)

# Plot activated terms
p1 <- emapplot(gseaGO_all_activated, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Activated Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Plot suppressed terms
p2 <- emapplot(gseaGO_all_suppressed, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Suppressed Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

p1 <- ridgeplot(gseaGO_all_simplified, showCategory = 20) + labs(x = "enrichment distribution")
print(p1)

gseaKEGG <- gseKEGG(geneList = genelist,
                organism = "hsa",
                keyType = "kegg",
                nPermSimple = 100000)
saveRDS(gseaKEGG, here::here(output_dir_data, "CM_GSEA_KEGG_mPAP.rds"))

gseaKEGG <- readRDS(here::here(output_dir_data, "CM_GSEA_KEGG_mPAP.rds"))
p1 <- dotplot(gseaKEGG, showCategory = 10, split=".sign", title = paste0("GSEA\nKEGG Pathways", "\n mPAP-associated genes"), label_format = 30, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"),
        strip.text = element_text(face = "bold", size = 18)) +
  facet_grid(.~.sign)
print(p1)

gseaKEGG@result <- classify_terms(gseaKEGG@result)
# Split the GSEA results
activated_terms <- subset(gseaKEGG@result, sign == "activated")
suppressed_terms <- subset(gseaKEGG@result, sign == "suppressed")

# Create two separate enrichResult objects
gseaKEGG_activated <- gseaKEGG
gseaKEGG_activated@result <- activated_terms

gseaKEGG_suppressed <- gseaKEGG
gseaKEGG_suppressed@result <- suppressed_terms

# Plot activated terms
p1 <- emapplot(pairwise_termsim(gseaKEGG_activated), showCategory = 20, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Activated Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Plot suppressed terms
p2 <- emapplot(pairwise_termsim(gseaKEGG_suppressed), showCategory = 20, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Suppressed Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

p1 <- ridgeplot(gseaKEGG, showCategory = 20) + labs(x = "enrichment distribution")
print(p1)

gseaGO_all_activated_df <- setReadable(gseaGO_all_activated, OrgDb = org.Hs.eg.db)
gseaGO_all_activated_df <- gseaGO_all_activated_df@result

Heatmaps

metadata_mPAP <- metadata[complete.cases(metadata[["mPAP_mmHg"]]),]

pseudocounts_mPAP <- pseudocounts[, metadata_mPAP$Sample]
pseudocounts_mPAP <- normLibSizes(pseudocounts_mPAP)

counts_mPAP <- cpm(pseudocounts_mPAP, log = TRUE)

## Set a color-blind friendly palette
heat_colors <- rev(brewer.pal(11, "PuOr"))

metadata_heatmap <- metadata_mPAP %>% dplyr::select(Sample, mPAP_mmHg, Batch_processing, Age_years) %>% dplyr::arrange(mPAP_mmHg)

DEGs_mPAP <- signif[["mPAP_mmHg"]]

sig_counts <- counts_mPAP[rownames(counts_mPAP) %in% DEGs_mPAP$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         show_rownames = FALSE,
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 10, 
         height = 20) 

Past versions of Heatmap_DEGs_mPAP-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_heatmap, sig_counts)

metadata_heatmap <- metadata_mPAP %>% dplyr::select(Sample, mPAP_mmHg, Batch_processing, Age_years) %>% dplyr::arrange(mPAP_mmHg)

DEGs_mPAP_top <- DEGs_mPAP
DEGs_mPAP_top <- DEGs_mPAP_top %>% 
  dplyr::filter(percentCells > 5) 
DEGs_mPAP_top_up <- DEGs_mPAP_top %>% 
  dplyr::filter(logFC > 2) %>% 
  dplyr::arrange(FDR)
DEGs_mPAP_top_down <- DEGs_mPAP_top %>% 
  dplyr::filter(logFC < -2) %>% 
  dplyr::arrange(FDR)

DEGs_mPAP_top <- dplyr::bind_rows(head(DEGs_mPAP_top_up, 20), head(DEGs_mPAP_top_down, 20))
rm(DEGs_mPAP_top_up, DEGs_mPAP_top_down)

sig_counts <- counts_mPAP[rownames(counts_mPAP) %in% DEGs_mPAP_top$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 15, 
         height = 20,
         show_rownames = T) 

Past versions of Heatmap_top_DEGs_mPAP-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_heatmap, sig_counts)

metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)

pseudocounts_mPAP_HC <- pseudocounts[, metadata_heatmap$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)

counts_mPAP_HC <- cpm(pseudocounts_mPAP_HC, log = TRUE)

sig_counts <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         show_rownames = FALSE,
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 10, 
         height = 20) 

Past versions of Heatmap_DEGs_HC-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

rm(metadata_heatmap, sig_counts)

metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)

pseudocounts_mPAP_HC <- pseudocounts[, metadata_heatmap$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)

counts_mPAP_HC <- cpm(pseudocounts_mPAP_HC, log = TRUE)

sig_counts <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP_top$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 15, 
         height = 20,
         show_rownames = T) 

Past versions of Heatmap_top_DEGs_HC-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

DEGs hierarchical clustering

# Filter and arrange metadata
metadata_mPAP_HC <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC") %>%
  dplyr::mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  dplyr::arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)

# Normalize pseudocounts and calculate log CPM
pseudocounts_mPAP_HC <- pseudocounts[, metadata_mPAP_HC$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)
counts_mPAP_HC <- edgeR::cpm(pseudocounts_mPAP_HC, log = TRUE)

# Filter counts matrix by DEGs
counts_matrix <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP$gene, ]

# Convert to long format and add Condition dividing DCM samples in 3 groups based on mPAP
counts_tibble <- counts_matrix %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "gene") %>%
  tidyr::pivot_longer(-gene, names_to = "sample", values_to = "log_pseudocounts") %>%
  dplyr::mutate(Condition = dplyr::case_when(
    sample %in% c("CZD4", "CZD5", "CZD6", "CZD7") ~ "HC",
    sample %in% c("PK439", "PK375", "PK229", "P5", "P8", "PK153") ~ "DCM - Mild PH",
    sample %in% c("PK385", "PK437", "PK323", "PK403", "PK213", "PK461") ~ "DCM - Moderate PH",
    sample %in% c("PK369", "PK357", "PK337", "PK345", "PK401", "PK463") ~ "DCM - Severe PH"
  ))

# Compute the mean log counts for each gene and condition
mean_tibble <- counts_tibble %>%
  dplyr::group_by(gene, Condition) %>%
  dplyr::summarise(mean_log_counts = mean(log_pseudocounts, na.rm = TRUE), .groups = 'drop')

# Convert to wide format
wide_mean_tibble <- mean_tibble %>%
  tidyr::pivot_wider(names_from = Condition, values_from = mean_log_counts)

# Convert the wide tibble to a matrix and transpose
mean_matrix <- wide_mean_tibble %>%
  tibble::column_to_rownames(var = "gene") %>%
  as.matrix()

mean_matrix  <- mean_matrix %>%
  t() %>% 
  scale() %>% 
  t()

# Perform hierarchical clustering
gene_dist <- dist(mean_matrix)
gene_hclust <- hclust(gene_dist, method = "complete")

# Plot the dendrogram
plot(gene_hclust, labels = FALSE)
abline(h = 3, col = "brown", lwd = 2) # Add horizontal line to illustrate cutting dendrogram

Past versions of Dendrogram-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

# Cut the dendrogram into 4 clusters
gene_cluster <- cutree(gene_hclust, k = 4) %>%
  tibble::enframe() %>%
  dplyr::rename(gene = name, cluster = value)

# Convert scaled matrix to long format and join with gene_cluster
mean_tibble_cluster <- mean_matrix %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "gene") %>%
  tidyr::pivot_longer(-gene, names_to = "Condition", values_to = "scaled_log_pseudocounts") %>%
  dplyr::inner_join(gene_cluster, by = "gene")

# Set factor levels for Condition
mean_tibble_cluster$Condition <- factor(mean_tibble_cluster$Condition, levels = c("HC", "DCM - Mild PH", "DCM - Moderate PH", "DCM - Severe PH"))

# Calculate the number of genes per cluster
cluster_counts <- mean_tibble_cluster %>%
  dplyr::group_by(cluster) %>%
  dplyr::summarise(gene_count = dplyr::n_distinct(gene))

# Create a named vector for custom facet labels
custom_labels <- setNames(paste("Cluster", cluster_counts$cluster, "-", cluster_counts$gene_count, "genes"), cluster_counts$cluster)

# Plot the data with custom facet labels
mean_tibble_cluster %>%
  ggplot(aes(Condition, scaled_log_pseudocounts)) +
  geom_line(aes(group = gene), alpha = 0.2) +
  geom_line(stat = "summary", fun = "median", colour = "brown", size = 1.5, aes(group = 1)) +
  facet_grid(cols = vars(cluster), labeller = as_labeller(custom_labels)) +
  theme(
    strip.text = element_text(face = "bold", size = 12),       # Make facet labels bold
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10) # Rotate x-axis labels by 45 degrees
  )

Past versions of Clustering_DEGs-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

Over-representation DEGs clusters

DEGs_clusters <- list()
DEGs_clusters[["cluster1"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "1")$gene)
DEGs_clusters[["cluster2"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "2")$gene)
DEGs_clusters[["cluster3"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "3")$gene)
DEGs_clusters[["cluster4"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "4")$gene)

DEGs_clusters_entrez <- list()
DEGs_clusters_entrez[["cluster1"]] <- bitr(DEGs_clusters[["cluster1"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster2"]] <- bitr(DEGs_clusters[["cluster2"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster3"]] <- bitr(DEGs_clusters[["cluster3"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster4"]] <- bitr(DEGs_clusters[["cluster4"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID

cluster1_GO <- enrichGO(DEGs_clusters_entrez[["cluster1"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster1_GO <- gsfilter(cluster1_GO, by = 'Count', min = 3)
nrow(cluster1_GO)

cluster2_GO <- enrichGO(DEGs_clusters_entrez[["cluster2"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster2_GO <- gsfilter(cluster2_GO, by = 'Count', min = 3)
nrow(cluster2_GO)

cluster3_GO <- enrichGO(DEGs_clusters_entrez[["cluster3"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster3_GO <- gsfilter(cluster3_GO, by = 'Count', min = 10)
nrow(cluster3_GO)
cluster3_GO <- simplify(
  cluster3_GO,
  cutoff = 0.7,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
cluster3_GO <- reorder_GO_by_pvalue(cluster3_GO)
saveRDS(cluster3_GO, here::here(output_dir_data, "CM_cluster3_ORA.rds"))

cluster3_GO <- readRDS(here::here(output_dir_data, "CM_cluster3_ORA.rds"))
dotplot(cluster3_GO, 
        showCategory = 15, 
        title = paste0("GO - Over-representation analysis", "\nDEGs cluster 3"), 
        label_format = 27, 
        font.size = 15) +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))

cluster4_GO <- enrichGO(DEGs_clusters_entrez[["cluster4"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster4_GO <- gsfilter(cluster4_GO, by = 'Count', min = 10)
nrow(cluster4_GO)
cluster4_GO <- simplify(
  cluster4_GO,
  cutoff = 0.7,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
cluster4_GO <- reorder_GO_by_pvalue(cluster4_GO)
saveRDS(cluster4_GO, here::here(output_dir_data, "CM_cluster4_ORA.rds"))

cluster4_GO <- readRDS(here::here(output_dir_data, "CM_cluster4_ORA.rds"))
dotplot(cluster4_GO, 
        showCategory = 15, 
        title = paste0("GO - Over-representation analysis", "\nDEGs cluster 4"), 
        label_format = 27, 
        font.size = 15) +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))

cluster3_GO <- pairwise_termsim(cluster3_GO)
cluster4_GO <- pairwise_termsim(cluster4_GO)

p1 <- emapplot(cluster3_GO, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("DEGs cluster 3") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Plot suppressed terms
p2 <- emapplot(cluster4_GO, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("DEGs cluster 4") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

Top genes

DEGs_mPAP_cluster1 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster1"]])
top_DEGs_mPAP_cluster1 <- DEGs_mPAP_cluster1 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(logFC) %>%
  head(20)

DEGs_mPAP_cluster2 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster2"]])
top_DEGs_mPAP_cluster2 <- DEGs_mPAP_cluster2 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(logFC) %>%
  head(20)

DEGs_mPAP_cluster3 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster3"]])
top_DEGs_mPAP_cluster3 <- DEGs_mPAP_cluster3 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(desc(logFC)) %>%
  head(20)

DEGs_mPAP_cluster4 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster4"]])
top_DEGs_mPAP_cluster4 <- DEGs_mPAP_cluster4 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(desc(logFC)) %>%
  head(20)

metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster1$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 1 DEGs") 

Past versions of Heatmap_DEGs_mPAP_cluster_1-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

ezInteractiveTableRmd(top_DEGs_mPAP_cluster1)

metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster2$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 2 DEGs") 

Past versions of Heatmap_DEGs_mPAP_cluster_2-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

ezInteractiveTableRmd(top_DEGs_mPAP_cluster2)

metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster3$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 3 DEGs") 

Past versions of Heatmap_DEGs_mPAP_cluster_3-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

ezInteractiveTableRmd(top_DEGs_mPAP_cluster3)

metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster4$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 4 DEGs") 

Past versions of Heatmap_DEGs_mPAP_cluster_4-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

ezInteractiveTableRmd(top_DEGs_mPAP_cluster4)

signif_cluster <- signif[["mPAP_mmHg"]] %>%
  left_join(gene_cluster) %>%
  left_join(ranked_genes)
write.csv(signif_cluster, file = here::here(output_dir_data, "CM_signif_cluster.csv"), quote=F, row.names = F)
ezInteractiveTableRmd(signif_cluster)

Scores

terms <- list(
  muscle_structure_development = c("actin", "muscle", "sarcomere", "cardiac", "fibril", "fiber", "actomyosin"),
  autophagy = c("autophag", "disassembly", "catabolic", "starvation"),
  transcription = c("transcription", "receptor"),
  vesicular_transport = c("transport", "Golgi", "vacuo", "vesicle", "endosom")
)

get_unique_genes <- function(df, terms) {
  filtered_df <- df %>%
    filter(grepl(paste(terms, collapse = "|"), Description, ignore.case = TRUE))
  unique(unlist(strsplit(filtered_df$core_enrichment, split = "/")))
}

genes_of_interest <- lapply(terms, get_unique_genes, df = gseaGO_all_activated_df)
names(genes_of_interest) <- names(terms)
saveRDS(genes_of_interest, here::here(output_dir_data, "CM_genes_of_interest.rds"))

for (i in 1:length(genes_of_interest)) {
  score_genes <- list(genes_of_interest[[i]])
  score_name <- paste0(names(genes_of_interest)[i], "_score")
  CM <- AddModuleScore(object = CM, features = score_genes, name = score_name)
}
metadata_CM <- CM@meta.data
colnames(metadata_CM) <- gsub("_score1", "_score", colnames(metadata_CM))
CM@meta.data <- metadata_CM
rm(metadata_CM)

plot_list <- lapply(paste0(names(genes_of_interest), "_score"), function(feature) {
  p <- FeaturePlot(CM, features = feature, order = FALSE)
  p + scale_color_viridis(discrete = FALSE, option = "turbo")
})

CombinePlots(plots = plot_list, ncol = 4)

Past versions of FeaturePlot_scores_1-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

plot_list <- lapply(paste0(names(genes_of_interest), "_score"), function(feature) {
  p <- FeaturePlot(CM, features = feature, order = TRUE)
  p + scale_color_viridis(discrete = FALSE, option = "turbo")
})

CombinePlots(plots = plot_list, ncol = 4)

Past versions of FeaturePlot_scores_2-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

VlnPlot(CM, features = paste0(names(genes_of_interest), "_score"), ncol = 4, group.by = "cell_state", pt.size = 0)

Past versions of VlnPlot_scores_cell_state-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

CM@meta.data[CM@meta.data$Sample %in% c("CZD4", "CZD5", "CZD6", "CZD7"), "Condition"] = "HC"
CM@meta.data[CM@meta.data$Sample %in% c("PK439", "PK375", "PK229", "P5", "P8", "PK153"), "Condition"] = "DCM - Mild PH"
CM@meta.data[CM@meta.data$Sample %in% c("PK385", "PK437", "PK323", "PK403", "PK213", "PK461"), "Condition"] = "DCM - Moderate PH"
CM@meta.data[CM@meta.data$Sample %in% c("PK369", "PK357", "PK337", "PK345", "PK401", "PK463"), "Condition"] = "DCM - Severe PH"
CM@meta.data[CM@meta.data$Sample %in% c("P3", "PK415", "PK441"), "Condition"] = "DCM - no mPAP data"

CM$Condition <- factor(CM$Condition, levels = c("HC", "DCM - Mild PH", "DCM - Moderate PH", "DCM - Severe PH", "DCM - no mPAP data"))

VlnPlot(CM, features = paste0(names(genes_of_interest), "_score"), ncol = 4, group.by = "Condition", pt.size = 0)

Past versions of VlnPlot_scores_Condition-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

Transcription factor inference analysis

net <- get_collectri(organism='human', split_complexes=FALSE)
net

# A tibble: 43,178 × 3
   source target   mor
   <chr>  <chr>  <dbl>
 1 MYC    TERT       1
 2 SPI1   BGLAP      1
 3 SMAD3  JUN        1
 4 SMAD4  JUN        1
 5 STAT5A IL2        1
 6 STAT5B IL2        1
 7 RELA   FAS        1
 8 WT1    NR0B1      1
 9 NR0B2  CASP1      1
10 SP1    ALDOA      1
# ℹ 43,168 more rows

ranked_genes_TF <- ranked_genes %>%
  dplyr::select(gene, logFC, metric, PValue) %>% 
  remove_rownames() %>%
  column_to_rownames(var = "gene") %>%
  as.matrix()
head(ranked_genes_TF)

               logFC      metric       PValue
A1CF      -4.4846666 -20.7823457 2.322259e-05
A2M       -1.6943794  -6.0929917 2.535118e-04
A2M-AS1    0.2481086   0.1532492 2.411737e-01
A2ML1     -0.4535097  -0.1800881 4.007756e-01
A2ML1-AS1 -2.5607983  -8.7931153 3.683495e-04
A4GALT     2.0319539   6.1047409 9.899886e-04

# Run ulm
contrast_acts <- run_ulm(mat=ranked_genes_TF[, 'metric', drop=FALSE], net=net, .source='source', .target='target',
                  .mor='mor', minsize = 5)
head(contrast_acts)

# A tibble: 6 × 5
  statistic source condition  score p_value
  <chr>     <chr>  <chr>      <dbl>   <dbl>
1 ulm       ABL1   metric     1.49    0.137
2 ulm       AHR    metric     1.06    0.287
3 ulm       AIP    metric     0.893   0.372
4 ulm       AIRE   metric    -0.730   0.465
5 ulm       AP1    metric     0.776   0.438
6 ulm       APEX1  metric     0.458   0.647

n_tfs <- 40

# Filter top TFs in both signs
f_contrast_acts <- contrast_acts %>%
  mutate(rnk = NA)
msk <- f_contrast_acts$score > 0
f_contrast_acts[msk, 'rnk'] <- rank(-f_contrast_acts[msk, 'score'])
f_contrast_acts[!msk, 'rnk'] <- rank(-abs(f_contrast_acts[!msk, 'score']))
tfs <- f_contrast_acts %>%
  arrange(rnk) %>%
  head(n_tfs) %>%
  pull(source)
f_contrast_acts <- f_contrast_acts %>%
  filter(source %in% tfs)

# Plot
ggplot(f_contrast_acts, aes(x = reorder(source, score), y = score)) + 
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", 
        mid = "whitesmoke", midpoint = 0) + 
    theme_minimal() +
    theme(axis.title = element_text(face = "bold", size = 12),
        axis.text.x = 
            element_text(angle = 45, hjust = 1, size =10, face= "bold"),
        axis.text.y = element_text(size =10, face= "bold"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
    xlab("TFs")

Past versions of mPAP_TF_barplot-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

FeaturePlot(CM, features = f_contrast_acts$source, ncol = 8)

Past versions of FeaturePlot_TF-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

contrast_acts_filtered <- contrast_acts %>%
  filter(source %in% filter(percent_stats, percentCells > 5)$gene)

n_tfs <- 40

# Filter top TFs in both signs
f_contrast_acts <- contrast_acts_filtered %>%
  mutate(rnk = NA) %>%
  filter(p_value < 0.05)
msk <- f_contrast_acts$score > 0
f_contrast_acts[msk, 'rnk'] <- rank(-f_contrast_acts[msk, 'score'])
f_contrast_acts[!msk, 'rnk'] <- rank(-abs(f_contrast_acts[!msk, 'score']))
tfs <- f_contrast_acts %>%
  arrange(rnk) %>%
  head(n_tfs) %>%
  pull(source)
f_contrast_acts <- f_contrast_acts %>%
  filter(source %in% tfs) 

# Plot
ggplot(f_contrast_acts, aes(x = reorder(source, score), y = score)) + 
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", 
        mid = "whitesmoke", midpoint = 0) + 
    theme_minimal() +
    theme(axis.title = element_text(face = "bold", size = 12),
        axis.text.x = 
            element_text(angle = 45, hjust = 1, size =10, face= "bold"),
        axis.text.y = element_text(size =10, face= "bold"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
    xlab("TFs")

Past versions of mPAP_TF_barplot_filtered-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

FeaturePlot(CM, features = f_contrast_acts$source, ncol = 6)

Past versions of FeaturePlot_TF_filtered-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

tf <- 'MYOCD'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Past versions of TF_MYOCD-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

tf <- 'SRF'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

tf <- 'MEF2A'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Past versions of TF_MEF2A-1.png

Version	Author	Date
182f9de	GinoBonazza	2024-08-19
5cd0f74	GinoBonazza	2024-08-01
61e8d17	GinoBonazza	2024-08-01

https://www.nature.com/articles/s41419-023-05665-8

tf <- 'MEF2C'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Past versions of TF_MEF2C-1.png

Version	Author	Date
182f9de	GinoBonazza	2024-08-19
61e8d17	GinoBonazza	2024-08-01

tf <- 'NFATC3'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Past versions of TF_NFATC3-1.png

Version	Author	Date
182f9de	GinoBonazza	2024-08-19
61e8d17	GinoBonazza	2024-08-01

# Run ulm
sample_acts <- run_ulm(mat=counts_mPAP_HC, net=net, .source='source', .target='target',
                  .mor='mor', minsize = 5)

# Transform to wide matrix
sample_acts_mat <- sample_acts %>%
  pivot_wider(id_cols = 'condition', names_from = 'source',
              values_from = 'score') %>%
  column_to_rownames('condition') %>%
  as.matrix()

# Get top tfs with more variable means across clusters
tfs <- arrange(f_contrast_acts, desc(score))$source
sample_acts_mat <- sample_acts_mat[,tfs]

# Scale per sample
sample_acts_mat <- scale(sample_acts_mat)
sample_acts_mat <- t(sample_acts_mat)
sample_acts_mat <- sample_acts_mat[, metadata_mPAP_HC$Sample, drop = FALSE]


# Choose color palette
palette_length = 100
my_color = colorRampPalette(c("Darkblue", "white","red"))(palette_length)

my_breaks <- c(seq(-3, 0, length.out=ceiling(palette_length/2) + 1),
               seq(0.05, 3, length.out=floor(palette_length/2)))

# Plot
pheatmap::pheatmap(sample_acts_mat, 
         border_color = NA, 
         color=my_color, 
         breaks = my_breaks,
         cluster_rows = F,
         cluster_cols = F,
         annotation_col = dplyr::select(metadata_mPAP_HC, -Batch_processing)) 

Past versions of TF_heatmap_Sample-1.png

Version	Author	Date
61e8d17	GinoBonazza	2024-08-01

saveRDS(CM, 
        here::here(output_dir_data, "CM.rds"))

Session information

sessionInfo()

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
 [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
 [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
[11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] htmltools_0.5.8.1           DT_0.32                    
 [3] ComplexHeatmap_2.18.0       UpSetR_1.4.0               
 [5] ezRun_3.19.1                Biostrings_2.70.2          
 [7] XVector_0.42.0              enrichplot_1.22.0          
 [9] dorothea_1.14.1             OmnipathR_3.10.1           
[11] decoupleR_2.9.7             EnhancedVolcano_1.20.0     
[13] ggrepel_0.9.5               speckle_1.2.0              
[15] edgeR_4.0.16                limma_3.58.1               
[17] ReactomePA_1.46.0           AnnotationHub_3.10.0       
[19] BiocFileCache_2.10.1        dbplyr_2.4.0               
[21] org.Hs.eg.db_3.18.0         AnnotationDbi_1.64.1       
[23] clusterProfiler_4.10.1      gprofiler2_0.2.3           
[25] harmony_1.2.0               clustree_0.5.1             
[27] ggraph_2.2.1                rhdf5_2.46.1               
[29] Rcpp_1.0.12                 data.table_1.15.2          
[31] plyr_1.8.9                  gtable_0.3.4               
[33] gtools_3.9.5                ArchR_1.0.2                
[35] statmod_1.5.0               patchwork_1.2.0            
[37] scater_1.30.1               scuttle_1.12.0             
[39] Matrix.utils_0.9.7          RColorBrewer_1.1-3         
[41] scales_1.3.0                knitr_1.45                 
[43] gridExtra_2.3               png_0.1-8                  
[45] pheatmap_1.0.12             reshape2_1.4.4             
[47] lubridate_1.9.3             forcats_1.0.0              
[49] stringr_1.5.1               purrr_1.0.2                
[51] tidyr_1.3.1                 tibble_3.2.1               
[53] tidyverse_2.0.0             magrittr_2.0.3             
[55] ggplot2_3.5.0               dplyr_1.1.4                
[57] scCustomize_2.1.2           scDblFinder_1.16.0         
[59] Matrix_1.6-5                DropletUtils_1.22.0        
[61] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[63] Biobase_2.62.0              GenomicRanges_1.54.1       
[65] GenomeInfoDb_1.38.7         IRanges_2.36.0             
[67] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[69] MatrixGenerics_1.14.0       matrixStats_1.2.0          
[71] SeuratObject_5.0.2          Seurat_4.4.0               
[73] xlsx_0.6.5                  readxl_1.4.3               
[75] readr_2.1.5                 here_1.0.1                 

loaded via a namespace (and not attached):
  [1] igraph_2.0.3                  graph_1.80.0                 
  [3] ica_1.0-3                     plotly_4.10.4                
  [5] rematch2_2.1.2                zlibbioc_1.48.0              
  [7] tidyselect_1.2.1              rvest_1.0.4                  
  [9] bit_4.0.5                     doParallel_1.0.17            
 [11] clue_0.3-65                   lattice_0.22-5               
 [13] rjson_0.2.21                  blob_1.2.4                   
 [15] S4Arrays_1.2.1                parallel_4.3.1               
 [17] cli_3.6.2                     ggplotify_0.1.2              
 [19] goftest_1.2-3                 BiocIO_1.12.0                
 [21] bluster_1.12.0                grr_0.9.5                    
 [23] BiocNeighbors_1.20.2          uwot_0.1.16                  
 [25] shadowtext_0.1.3              curl_5.2.1                   
 [27] mime_0.12                     evaluate_0.23                
 [29] tidytree_0.4.6                leiden_0.4.3.1               
 [31] stringi_1.8.3                 backports_1.4.1              
 [33] XML_3.99-0.16.1               httpuv_1.6.14                
 [35] paletteer_1.6.0               rappdirs_0.3.3               
 [37] splines_4.3.1                 logger_0.3.0                 
 [39] bcellViper_1.38.0             sctransform_0.4.1            
 [41] ggbeeswarm_0.7.2              DBI_1.2.2                    
 [43] HDF5Array_1.30.1              jquerylib_0.1.4              
 [45] reactome.db_1.86.2            withr_3.0.0                  
 [47] git2r_0.33.0                  rprojroot_2.0.4              
 [49] xgboost_1.7.7.1               lmtest_0.9-40                
 [51] ggnewscale_0.4.10             tidygraph_1.3.1              
 [53] rtracklayer_1.62.0            BiocManager_1.30.22          
 [55] htmlwidgets_1.6.4             fs_1.6.3                     
 [57] labeling_0.4.3                SparseArray_1.2.4            
 [59] cellranger_1.1.0              reticulate_1.35.0            
 [61] zoo_1.8-12                    timechange_0.3.0             
 [63] foreach_1.5.2                 fansi_1.0.6                  
 [65] ggtree_3.10.1                 R.oo_1.26.0                  
 [67] irlba_2.3.5.1                 ggrastr_1.0.2                
 [69] gridGraphics_0.5-1            ellipsis_0.3.2               
 [71] lazyeval_0.2.2                yaml_2.3.8                   
 [73] survival_3.5-8                scattermore_1.2              
 [75] BiocVersion_3.18.1            crayon_1.5.2                 
 [77] RcppAnnoy_0.0.22              progressr_0.14.0             
 [79] tweenr_2.0.3                  later_1.3.2                  
 [81] ggridges_0.5.6                codetools_0.2-19             
 [83] GlobalOptions_0.1.2           KEGGREST_1.42.0              
 [85] Rtsne_0.17                    shape_1.4.6.1                
 [87] Rsamtools_2.18.0              filelock_1.0.3               
 [89] pkgconfig_2.0.3               xml2_1.3.6                   
 [91] GenomicAlignments_1.38.2      aplot_0.2.2                  
 [93] spatstat.sparse_3.0-3         ape_5.7-1                    
 [95] viridisLite_0.4.2             xtable_1.8-4                 
 [97] highr_0.10                    httr_1.4.7                   
 [99] tools_4.3.1                   globals_0.16.3               
[101] beeswarm_0.4.0                checkmate_2.3.1              
[103] nlme_3.1-164                  selectr_0.4-2                
[105] HDO.db_0.99.1                 crosstalk_1.2.1              
[107] digest_0.6.35                 farver_2.1.1                 
[109] tzdb_0.4.0                    yulab.utils_0.1.4            
[111] viridis_0.6.5                 glue_1.7.0                   
[113] cachem_1.0.8                  polyclip_1.10-6              
[115] generics_0.1.3                parallelly_1.37.1            
[117] ScaledMatrix_1.10.0           pbapply_1.7-2                
[119] vroom_1.6.5                   spam_2.10-0                  
[121] gson_0.1.0                    dqrng_0.3.2                  
[123] utf8_1.2.4                    graphlayouts_1.1.1           
[125] shiny_1.8.0                   GenomeInfoDbData_1.2.11      
[127] R.utils_2.12.3                rhdf5filters_1.14.1          
[129] RCurl_1.98-1.14               memoise_2.0.1                
[131] rmarkdown_2.26                R.methodsS3_1.8.2            
[133] future_1.33.1                 RANN_2.6.1                   
[135] Cairo_1.6-2                   spatstat.data_3.0-4          
[137] rstudioapi_0.15.0             cluster_2.1.6                
[139] whisker_0.4.1                 janitor_2.2.0                
[141] spatstat.utils_3.0-4          hms_1.1.3                    
[143] fitdistrplus_1.1-11           munsell_0.5.0                
[145] cowplot_1.1.3                 colorspace_2.1-0             
[147] rlang_1.1.4                   DelayedMatrixStats_1.24.0    
[149] sparseMatrixStats_1.14.0      dotCall64_1.1-1              
[151] ggforce_0.4.2                 circlize_0.4.16              
[153] xfun_0.42                     iterators_1.0.14             
[155] abind_1.4-5                   GOSemSim_2.28.1              
[157] interactiveDisplayBase_1.40.0 treeio_1.26.0                
[159] rJava_1.0-11                  Rhdf5lib_1.24.2              
[161] bitops_1.0-7                  promises_1.2.1               
[163] scatterpie_0.2.1              RSQLite_2.3.5                
[165] qvalue_2.34.0                 fgsea_1.28.0                 
[167] DelayedArray_0.28.0           GO.db_3.18.0                 
[169] compiler_4.3.1                prettyunits_1.2.0            
[171] beachmat_2.18.1               graphite_1.48.0              
[173] listenv_0.9.1                 workflowr_1.7.1              
[175] BiocSingular_1.18.0           tensor_1.5                   
[177] MASS_7.3-60.0.1               progress_1.2.3               
[179] BiocParallel_1.36.0           spatstat.random_3.2-3        
[181] R6_2.5.1                      fastmap_1.1.1                
[183] fastmatch_1.1-4               vipor_0.4.7                  
[185] ROCR_1.0-11                   rsvd_1.0.5                   
[187] KernSmooth_2.23-22            miniUI_0.1.1.1               
[189] deldir_2.0-4                  bit64_4.0.5                  
[191] spatstat.explore_3.2-6        lifecycle_1.0.4              
[193] ggprism_1.0.4                 restfulr_0.0.15              
[195] xlsxjars_0.6.1                sass_0.4.9                   
[197] vctrs_0.6.5                   spatstat.geom_3.2-9          
[199] snakecase_0.11.1              DOSE_3.28.2                  
[201] scran_1.30.2                  ggfun_0.1.4                  
[203] sp_2.1-3                      future.apply_1.11.1          
[205] bslib_0.6.1                   pillar_1.9.0                 
[207] metapod_1.10.1                locfit_1.5-9.9               
[209] jsonlite_1.8.8                GetoptLong_1.0.5